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Abstract 

We  discuss  domain  decomposition  methods  with  which  the  often 
very  large  hnear  systems  of  algebraic  equations,  arising  when  elliptic 
problems  are  discretized  by  finite  differences  or  finite  elements,  can 
be  solved  with  the  aid  of  exact  or  approximate  solvers  for  the  same 
equations  restricted  to  subregions.  The  interaction  between  the  subre- 
gions,  to  enforce  appropriate  continuity  requirements,  is  handled  by  an 
iterative  method,  often  a  preconditioned  conjugate  gradient  method. 
Much  of  the  work  is  local  and  cjin  be  carried  out  in  parallel.  We  first 
explore  how  ideas  from  structural  engineering  computations  naturally 
lead  to  certain  matrix  splittings.  In  preparation  for  the  detailed  design 
and  analysis  of  related  domain  decomposition  methods,  we  then  con- 
sider the  Schwarz  alternating  algorithm,  discovered  already  in  1869. 
That  algorithm  can  conveniently  be  expressed  in  terms  of  certain  pro- 
jections. We  develop  these  ideas  further  and  discuss  an  interesting 
additive  variant  of  the  Schwarz  method.  This  also  leads  to  the  devel- 
opment of  a  general  framework,  which  already  has  proven  quite  useful 
in  the  study  of  a  variety  of  domain  decomposition  methods  and  certain 
related  algorithms.  We  demonstrate  this  by  developing  several  algo- 
rithms and  by  showing  how  their  rates  of  convergence  can  be  estimated. 
One  of  them  is  a  Schwarz- type  method,  for  which  the  subregions  over- 
lap, while  the  others  aje  so  called  iterative  substructuring  methods, 
where  the  subregions  do  not  overlap.  Compared  to  previous  studies 
of  iterative  substructuring  methods,  our  proof  is  simpler  and  in  one 
case  it  can  be  completed  without  using  a  finite  element  extension  the- 
orem. Such  a  theorem  has,  to  our  knowledge,  always  been  used  in  the 
previous  ameJysis  in  all  but  the  very  simplest  cases. 


1      Introduction. 

Domain  decomposition  methods  have  recently  become  an  important  focus 
of  research  on  numerical  methods  for  partial  differential  equations.  They 
appear  to  offer  the  best  promise  for  the  parallel  solution  of  the  often  very 
large  systems  of  linear  or  nonlinear  algebraic  systems  of  equations  which 
arise  when  the  elliptic  problems  of  elasticity,  fluid  dynamics  and  many  other 
important  applications  axe  discretized  by  finite  elements  or  finite  differences. 
The  rapidly  growing  interest  in  this  field  is  reflected  in  a  series  of  SIAM 
sponsored  symposia  and  mini-symposia  exclusively  devoted  to  research  in 
this  area;  see  [9,15]. 

The  domain  decomposition  methods  considered  here  can  be  regarded  as 
divide  and  conquer  cilgorithms.  In  each  step  of  an  iteration,  the  original 
discrete  elliptic  problem  is  solved  on  the  subregions  into  which  the  origi- 
nal region  has  been  divided.  (Under  favorable  circumstances,  an  already 
existing  code  can  be  used  for  this  purpose.)  The  interaction  between  the 
different  parts  of  the  region  is  handled  by  transferring  suitable  data  across 
the  interfaces,  which  are  created  by  the  subdivision.  These  data  are  gener- 
ated in  each  iteration  from  the  residual  of  the  original  or  a  derived  system 
of  equations.  A  conjugate  gradient  method  is  often  used  to  accelerate  the 
convergence.  In  addition  to  the  local  problems,  some  global,  coarse  problem 
must  be  incorporated  into  the  preconditioner  in  order  to  obtain  a  fast  rate 
of  convergence  in  the  case  of  many  subregions.  That  is  the  case  of  primary 
interest  in  parallel  computing  research.  We  note  that  the  need  for  a  coarse, 
global  part  of  the  preconditioner  is  well  known  in  multi-grid  work.  More  for- 
maUy,  we  have  shown  in  our  previous  work  that  if  in  each  iteration  step  of 
a  domain  decomposition  algorithm  information  is  ordy  exchanged  between 
neighboring  substructures,  the  rate  of  convergence  can  be  no  better  than 
if  the  conjugate  gradient  method  is  used,  without  any  preconditioning,  for 
the  coarse  model  obtained  by  using  the  subregions  as  elements;  cf.  Widlund 
[29]. 

There  has  been  considerable  progress  in  this  research  area  and  a  number 
of  fast  algorithms  have  been  designed  and  studied  for  which  the  condition 
number  of  the  iteration  matrix  is  uniformly  bounded  or  grows  only  in  pro- 
portion to  (1  -I-  log{n/h)y,q  =  2  or  3,  where  H  is  the  diameter  of  a  typical 
subregion  and  h  the  diameter  of  a  typical  element  into  which  the  subregions 
are  divided;  see  e.g.  Bramble,  Pasciak  and  Schatz  [5,6,7,8]  ,  Dryja  [12], 
Dryja,  Proskurowski  and  Widlund  [13],  Dryja  and  Widlund  [14]  and  Wid- 
lund [29]  .  Most  of  the  important  results  in  this  field  have  been  developed 


in  a  finite  element  framework  and  we  also  adopt  such  an  approach  in  this 
paper. 

To  assess  the  complexity  of  domain  decomposition  algorithms,  tradi- 
tional tools  must  also  be  used  to  measure  the  cost  of  solving  the  different 
subproblems.  If  a  parallel  computing  system  is  used,  there  are  of  course 
many  additional  aspects  to  estimating  the  overall  cost  of  the  computation. 
We  provide  no  detailed  discussion  of  such  matters  here.  Instead,  we  focus  on 
the  design  of  preconditioners  that  give  rapid  rates  of  convergence  while  de- 
creasing the  amount  of  work  per  step,  and  on  the  development  of  a  general 
framework  inside  which  a  variety  of  algorithms  can  be  designed  and  ana- 
lyzed. We  note  that  several  of  sets  of  careful  numerical  experiments  with 
various  domain  decomposition  algorithms  have  been  carried  out  on  parallel 
computers;  cf.  Keyes  and  Gropp  [18,19],  who  used  a  hypercube,  and  Green- 
baum  et  al.  [17],  who  used  an  ultra  computer  prototype  built  at  the  Courant 
Institute. 

Domain  decomposition  methods  can  be  classified  according  to  whether 
the  subregions  (substructures)  overlap  or  not.  The  algorithms  which  do 
not  use  any  overlap  are  often  called  iterative  substructuring  methods  while 
the  others  are  called  Schwarz-type  methods.  This  is  in  recognition  of  the 
importance  of  ideas  of  structural  engineering  computing  in  the  development 
of  the  former  class  of  methods  and  of  the  fundamental  contribution  by  H.  A. 
Schwarz,  who  introduced  his  alternating  method  in  1869;  cf.  Schwarz  [24]. 
The  two  classes  of  methods  have  much  in  common  and  there  is  a  real  promise 
that  a  unified  theory  can  be  developed.  This  paper  is  aa  effort  in  that 
direction.  We  are  able  to  show  that  some  interesting  iterative  substructuring 
methods  can  also  be  derived  from  a  so  called  additive  variant  of  the  Schwarz 
algorithm.  Technically  this  work  also  offers  something  new.  One  of  the 
algorithms  can  be  analyzed  without  using  a  so  called  finite  element  extension 
theorem.  Such  theorems  have  previously  been  central  to  the  development 
of  the  theory;  see  Widlund  [28]  for  a  proof  of  such  a  theorem  for  general 
conforming  finite  element  methods  and  a  general  discussion. 

We  note  that  in  a  recent  paper,  Bjorstad  and  Widlund  [3]  have  returned 
to  the  two  subregion  case  to  show  that  a  method,  introduced  by  Chan  and 
Resasco  [10],  is  in  fact  the  classical  Schwarz'  method  accelerated  by  the 
conjugate  gradient  method.  That  study  provides  additional  evidence  that 
we  need  not  strictly  distinguish  between  methods  which  use  overlapping 
subregions  and  those  who  do  not. 

In  our  development  of  the  theory,  a  prominent  role  is  played  by  sub- 
spaces  and  projections.  Already  fifty  years  ago  Sobolev  [25]  showed  that  the 


classical  Schwarz  algorithm  can  be  analyzed  using  a  variational  framework 
and  more  recently  this  approach  has  been  further  developed  by  Pierre- Louis 
Lions  [20].  An  additive  Schwarz  algorithms  was  first  studied  in  Dryja  and 
Widlund  [11].  In  this  paper,  we  strengthen  and  simplify  that  result.  We 
note  that  similar  ideas  have  also  been  discussed  by  Matsokin  and  Nepom- 
nyaschikh  [22].  The  general  ancdytic  framework  that  is  now  available  for  the 
study  of  these  additive  methods,  has  also  proven  very  useful  in  our  study  of 
so  called  iterative  refinement  methods;  cf.  Widlund  [32,33]. 

A  few  years  ago,  much  of  the  work  on  iterative  substructuring  methods 
was  focused  on  elliptic  problems  defined  on  regions  divided  into  two  or  a  few 
subregions  with  separating  curves  (surfaces)  which  do  not  intersect  ;  see  e.g. 
Bj0rstad  and  Widlund  [2]  and  Widlund  [28].  In  this  paper,  we  show  that 
some  of  those  ctlgorithms  and  results  can  be  recycled  and  combined  with 
an  iterative  substructuring  algorithm  derived  by  using  subspaces  and  pro- 
jections to  obtain  previously  known  algorithms  and  results.  This  approach 
highlights  how  preconditioners  can  be  built  from  parts  which  are  strictly 
local  to  cin  individual  substructure,  parts  which  involve  interaction  between 
pairs  of  neighboring  substructures  and  a  coarse  global  model  with  relatively 
few  degrees  of  freedom. 

This  paper  is  organized  as  follows.  In  Section  2,  we  review  some  of 
the  idecis  of  substructuring  that  are  very  important  in  the  development  of 
computational  methods  of  structural  engineering.  This  discussion  naturally 
leads  to  matrix  splittings,  which  provide  preconditioners  for  the  large  lin- 
ear systems  of  algebraic  equations,  which  arises  in  finite  element  work.  In 
section  3,  we  discuss  different  Schwarz  methods  and  some  general  tools  for 
estimating  their  rates  of  convergence.  In  the  concluding  sections,  we  show 
how  two  types  of  domain  decomposition  algorithms  can  be  analyzed  by  using 
relatively  simple  tools  of  mathematical  and  finite  element  analysis.  While 
we  can  do  much  with  linear  algebra,  we  ultimately  have  to  resort  to  tools  of 
ancilysis  in  order  to  complete  the  proofs  of  our  main  results. 


2      Substructures,  subspaces  and  projections. 

The  domain  decomposition  methods,  considered  in  this  paper,  provide  pre- 
conditioners of  systems  of  equations  representing  discrete  elliptic  problems. 
Preconditioners  can  be  associated  with  matrix  splittings  but  we  note  that 
those  considered  here  differ  in  certain  respects  from  most  of  those  encoun- 
tered the  standard  literature  on  iterative  methods;  cf  Varga  [27].  Domain  de- 


composition  methods  can  best  be  understood  in  terms  of  the  substructuring 
ideas,  which  have  provided  a  very  successful  framework  for  the  development 
of  very  large  programming  systems  for  structural  mechanics  computing;  cf. 
Bell,  Hatlestad,  Hansteen  and  Araldsen  [1]  or  Przemieniecki  [23]. 

We  consider  a  linear,  self  adjoint,  elliptic  problem,  which  is  discretized 
by  a  finite  element  method  on  a  bounded  Lipschitz  region.  The  region  il 
is  a  subset  of  iZ",  n=2  or  3,  the  differential  operator  is  the  Laplacian,  and 
zero  Dirichlet  conditions  and  continuous,  piecewise  linear  finite  elements  are 
used.  The  theory  could  equally  well  be  developed  for  much  more  general 
linear  elliptic  problems,  which  can  be  formulated  as  minimization  problems. 
Arbitrary  conforming  finite  elements  could  also  be  considered  without  fur- 
ther major  complications.  Nonconforming  finite  elements,  non-self  adjoint 
problems  and  problems  that  give  rise  to  indefinite  symmetric  systems  of 
equations  are  also  quite  important.  Some  progress  has  already  been  made 
in  such  cases. 

In  the  model  case  considered  here,  the  continuous  and  discrete  problems 
take  the  form 

a(u,r)  =  /(t;),Vt;e   V, 

and 

a{uH,Vh)  =  f{vH),^Vf,e    V\  (1) 

respectively.  The  bilinear  form  is  defined  by 

Vu  •  Vv  dx  . 


i(u,t;)=  / 
Jn 


This  form  defines  a  semi-norm  |«|//i(fi)  =  (a(")  w))*^^  in  the  Sobolev  space 
H\Q).  It  is  a  norm  for  V  =  n^{Q).  Here  H^i^)  is  the  subspace  of  H^^) 
functions  with  zero  trace  ;  all  elements  of  V  and  its  subspace  V  vanish 
on  dil,  the  boundary  of  il.  The  triangulation  of  fi  is  introduced  in  the 
following  way.  The  region  is  divided  into  non-overlapping  substructures 
n,  ,i  =  l,---,iV  .  To  simplify  the  description,  our  study  is  confined  to 
triangular  (simplicial)  substructures.  In  such  a  case,  the  original  region  must 
of  course  be  a  polygon  (polyhedron).  AH  the  substructures  f2,  are  further 
divided  into  elements.  The  common  assumption  in  finite  element  theory 
that  all  elements  are  shape  regular  is  adopted  and  the  sajne  assumption  is 
made  concerning  the  substructures.  On  the  element  level  this  means  that 
there  is  a  bound  on  HkI Pk  i  which  is  independent  of  the  number  of  degrees 
of  freedom  and  of  K .  Here  hj^  is  the  diaimeter  of  the  element  K  and  px  the 
diameter  of  the  largest  sphere  that  can  be  inscribed  in  K. 


Since  a{uh,Vh)  =  a(u,Vh),  V  v/,  G  V*  ,  the  finite  element  solution  is  the 
projection  of  the  exact  solution  onto  the  finite  element  space  with  respect  to 
the  inner  product  defined  by  the  bilinear  form.  We  will  see  that  the  problems 
defined  on  the  subregions,  from  which  preconditioners  for  the  entire  problem 
can  be  assembled,  can  similarly  be  viewed  in  terms  of  orthogonal  projections, 
onto  subspaces  directly  associated  with  the  subregion  in  question. 

By  using  ideas  of  structural  engineering  ,  the  stiffness  matrix  K  can  be 
constructed  in  the  following  way.  The  elements  of  K  are  given  by 

A:.,j  =  a(v.,V;)  , 

where  (p,  and  <^j  are  standard  finite  element  basis  functions.  Since  an  in- 
tegral over  n  can  be  written  as  a  sum  of  integrals  over  the  substructures, 
the  stiffness  matrix  can  be  assembled  from  the  stiffness  matrices  A'^'^'  the 
elements  of  which  are  defined  by 

with  i  and  j  corresponding  to  the  degrees  of  freedom  the  substructure  Q^. 
The  form  aQ^{uh,Vh)  represents  the  contribution  to  the  integral  from  that 
substructure. 

This  so  called  subassembly  process  can  be  summarized  in  the  formula 

where  i^')  is  the  subvector  of  parameter  values  associated  with  the  substruc- 
ture fi,  and  its  boundary  50,. 

In  engineering  computing  practice,  the  large  linear  system  with  the  co- 
efficient matrix  K  is  solved  using  programs,  which  are  elaborate  implemen- 
tations of  block  Gaussian  elimination;  cf.  e.g.  [1].  Particular  attention  is 
placed  on  the  use  of  data  structures  etc.,  which  allows  for  efficient  I/O,  and 
on  the  need  to  integrate  the  factorization  and  solution  steps  as  much  as 
possible  with  the  generation  of  the  blocks  which  together  define  the  stiffness 
matrix.  By  the  so-called  subassembly  process,  described  above  in  a  spe- 
cial case,  the  contributions  from  the  individual  elements  are  first  computed 
and  these  element  stiffness  matrices  are  then  merged  with  their  neighbors, 
creating  so-called  super  elements.  This  process  is  often  used  recursively  a 
substantial  number  of  times.  On  each  level,  the  variables  associated  with 
the  formation  of  a  particular  super  element  can  be  divided  into  two  sets; 
those  which  are  common  to  other  super  elements  and  the  interior  variables 


which  axe  not.  In  our  description  of  domain  decomposition  algorithms  only 
three  levels  are  considered,  the  elements,  with  a  characteristic  diameter  h, 
the  substructures  with  a  diameter  H  and  the  entire  region  fi,  which,  without 
loss  of  generality,  is  assumed  to  have  unit  diameter. 

The  last  phases  of  a  standard  engineering  calciilation  of  this  type,  involve 
a  substantial  fraction  of  the  arithmetical  work,  and  require  more  global 
communication  than  in  the  previous  phases,  where  the  work  can  proceed  in 
a  distributed  fashion,  without  synchronization  between  the  substructures. 
The  iterative  substructuring  methods  require  synchronization  once  every 
iteration,  but  the  best  of  them  converge  quite  fast.  They  also  demand 
much  less  communication  of  data  between  the  substructures  than  methods 
which  use  direct  methods  throughout.  The  use  of  inexact  solvers  for  the 
subproblems  can  also  be  considered  when  an  iterative  substructuring  method 
is  used.  If  a  substructure  has  a  simple  geometry,  special  fast  solvers  may 
also  be  used. 

If  we  divide  the  subvectors  x^'^  associated  with  the  i-th  substructure 
into  two,  Xj  and  Xg,  corresponding  to  the  variables  shared  with  other 
substructures  and  those  which  are  interior  to  the  substructure,  then  the 
matrix  A'^'^  can  be  written  as 

i  A-!;'   Kfs  \ 

Since  the  interior  variables  axe  associated  with  only  one  of  the  substructures, 
they  can  be  eliminated  locally  and  in  parallel.  The  reduced  matrix  is  a  so- 
called  Schur  complement  and  has  the  form 

^        -  ^BB-  ^IB    ^11       ^IB- 
It  is  now  easy  to  show  that  if  the  corresponding  Schur  complement  of  the 
global  stiffness  matrix  K  is  denoted  by  5,  then 

,TSy=^,fs^^h^^.  (2) 

The  elimination  of  the  interior  variables  from  the  substructures  can 
be  viewed  in  terms  of  orthogonal  projections,  with  respect  to  the  bilinear 
form,  of  the  solution  u^  of  equation  1  onto  the  subspaces  ^o(^<)  H  ^^i  '  = 
1,  •  •  • ,  JV.  These  subspaces  can  easily  be  shown  to  be  orthogonal,  in  the  sense 
of  the  bilinezo'  form,  to  the  so-called  piecewise  discrete  haurmonic  functions 
which  satisfy 

i^i?xP  +  A'gxg)  =  0,    Vi. 


If  the  locaJ  problems  are  solved  exactly  ,  what  remains  is  to  find  a  suffi- 
ciently accurate  approximation  of  the  part  of  the  solution  which  is  piecewise 
discrete  harmonic.  This  is  done  by  approximately  solving  the  reduced  linear 
system  with  the  matrix  S.  An  iterative  substructuring  method  is  obtained 
by  selecting  a  preconditioner  for  the  the  matrix  S.  Once  an  approximation 
of  the  solution  hcis  been  found  on  the  boundaries  of  the  substructures,  the 
solution  can  be  found  everywhere  by  solving  local  Dirichlet  problems  on  each 
substructure  separately. 

The  matrices  5'*'  are  full.  Complete  information  on  the  sparsity  and 
block  structure  of  S  can  be  obtained  by  using  equation  2.  It  is  also  natu- 
ral to  partition  the  vector  Xg  and  the  Schur  complement  S^'^  further.  We 
group  the  variables  associated  with  the  interior  of  each  edge  (face)  of  the 
substructure  into  separate  subvectors.  Each  of  them  is  associated  with  ex- 
actly two  substructures  which  are  neighbors  of  each  other.  The  remaining 
degrees  of  freedom  of  Xg  form  a  separate  subvector.  In  two  dimensions,  the 
components  of  this  vector,  Xy  ,  are  the  the  values  at  the  vertices  of  17,  ,  while 
in  three  dimensions  they  are  associated  with  all  the  nodes  of  the  edges  of  the 
simplex,  including  its  vertices.  We  can  say  that  these  are  the  nodes  of  the 
wire  basket  outlining  the  substructure.  The  Schur  complement  5'*'  is  thus 
represented  by  a  four  by  four  (five  by  five)  block  matrix.  If  one  or  several 
edges  of  il,  is  part  of  dil,  and  thus  not  associated  with  any  degrees  of  free- 
dom in  the  Dirichlet  case,  then  the  numbers  of  blocks  of  S^'^  is  reduced.  In 
what  follows,  we  consider  splittings  where  the  off-  diagonal  blocks  of  5^''  are 
dropped,  and  others  where  all  the  contributions  from  certain  substructures 
are  left  out  in  the  construction  of  the  preconditioner. 

In  the  case  of  two  substructures,  there  is  only  one  set  of  interface  vari- 
ables and  S  is  the  sum  of  S'^^  and  5'^^  the  Schur  complements  originating 
from  the  two  subregions  respectively.  The  basic  Neumann-Dirichlet  algo- 
rithm amounts  to  using  5'^)  or  5^^^  as  a  preconditioner  for  5;  of.  Bjorstad 
and  Widlund  [2]  and  Widlund  [28].  This  algorithm  can  be  extended  directly 
to  the  case  of  many  substructures  if  a  red-black  ordering  of  the  substructures 
exists;  cf.  Dryja,  Proskurowski  and  Widlund  [13]  and  Widlund  [29].  This 
splitting  corresponds  to  dropping  all  contributions  from  the  red  (or  black) 
substructures  in  the  sum  given  in  equation  2.  The  resulting  preconditioner 
corresponds  to  a  globed  problem  since  the  remaining  black  substructures  are 
connected  at  their  vertices  (wire  baskets).  If  we  group  all  edge  (face)  vari- 
ables and  all  the  vertex  (wire  basket)  variables  into  two  subvectors,  then  we 


can  write  5^*^  as 


( 


c(«)  c(«) 

^EE  "^EV 

Mi)T  Mi) 

■->EV  ■^VV 


In  the  quadratic  form,  which  represents  the  sum  over  the  red  substructures 
only,  the  variables  of  i^'  appear  only  once.  They  are  coupled  only  to  vari- 
ables which  are  associated  with  the  same  substructure.  The  system  of  linear 
equations,  which  correspond  to  the  Neumann-Dirichlet  preconditioner,  can 
therefore  be  written  in  terms  of  a  two  by  two  block  matrix  where  the  leading 
block  is  the  direct  sum  of  the  5^^;  which  correspond  to  the  red  substruc- 
tures.   An  analysis  of  the  structure  of  this  large  matrix  reveals  that  it  is 
quite  economical  to  use  standard  block  Gaussian  elimination,  in  the  case  of 
two  dimensions.  The  Schur  complement  that  remains  after  the  elimination 
of  the  edge  variables  is  sparse  and  its  dimension  is  equal  to  the  number  of 
variables  associated  with  substructure  vertices  which  do  not  fall  on  dO.  ;  see 
further  Dryja,  Proskurowski  and  Widlund  [13]  .    It  is  shown  in  Widlund 
[29],  that  the  condition  number  of  the  resulting  preconditioner  grows  in  pro- 
portion to  (1  -I-  log{n/h)y.    The  algorithm  is  less  attractive  in  the  three 
dimensional  case,  since  the  number  of  variables,  which  have  to  be  treated 
in  a  special  way,  are  relatively  more  numerous  than  for  a  problem  in  the 
plane.  However,  the  same  kind  of  bound  can  be  established  on  the  rate  of 
convergence.  There  is  an  interesting  variant  of  the  algorithm,  where  the  wire 
basket  variables  are  handled  differently  and  the  preconditioner  is  cheaper  to 
use;  see  Dryja  [12]  .  We  will  not  discuss  it  further  in  this  paper. 

Other  splittings  of  5  are  known,  which  give  as  good  (or  better)  results 
as  the  Neumann-Dirichlet  method.  We  note  that  when  that  algorithm  was 
considered,  the  values  at  the  vertices  (on  the  wire  basket)  were  treated  differ- 
ently from  those  of  the  interior  of  the  edges  (faces) .  The  methods  considered 
by  Bramble,  Pasciak  and  Schatz  [5,8]  also  treat  these  sets  of  variables  dif- 
ferently. As  was  shown  in  Widlund  [29,30],  the  part  of  the  preconditioner 
introduced  in  [5],  that  relates  to  the  vertex  variables  of  plane  problems,  can 
be  viewed  in  terms  of  a  projection  onto  a  finite  element  subspace,  where 
the  substructures  play  the  role  of  elements.  The  dimension  of  this  subspace, 
V^,  is  equal  to  the  number  of  substructure  vertices  which  belong  to  the 
open  set  Cl  .  In  the  splitting  which  corresponds  to  this  preconditioner,  the 
couplings  between  the  different  edges  (faces)  are  also  ignored.  Thus  in  the 
four  by  four  blocks  which  represent  the  5^'^,  the  off  diagonal  blocks  are  set 
to  zero.  The  diagonal  blocks  are  also  changed;  cf.  Section  5  .  In  Section  5, 
we  will  show  that  the  variables  associated  with  an  individual  edge  can  be 
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naturally  related  to  a  subspace  of  functions  which  vanish  outside  the  two 
substructures,  which  have  this  edge  in  common.  The  few  preconditioners 
which  have  been  studied  in  detail  for  the  three  dimensional  case,  cf.  Bram- 
ble, Pasciak  and  Schatz  [8]  ajid  Dryja  [12],  can  also  be  described  in  similar 
general  terms. 

3      Schwarz  Methods. 

We  begin  by  briefly  discussing  the  classical  formulation  of  Schwarz'  method 
in  the  continuous  case.  There  aje  two  fractional  steps  corresponding  to  two 
overlapping  subregions,  Q\  and  fiji  the  union  of  which  is  the  region  fi  .  The 
initial  guess  u°  G  V.  The  iterate  u"+^  is  determined  from  u"  by  sequentially 
updating  the  approximate  solution  on  the  two  subregions. 


aind 


„n+i     ^     „„+i/2  on  on'. 


We  could  equally  well  have  written  down  the  finite  element  version  of 
the  algorithm.  From  now  on,  we  only  consider  that  case.  It  is  easy  and 
convenient  to  describe  this  classical  method  in  terms  of  two  projections 
P„  i  =  1,2,  onto  v.'*  =  n^{n',)f]V''  ■  cf.  Lions  [20].  They  are  defined  by 

a{P,Vh,4>k)  =  a{vh,4>h),  ^4>H  €  V.''  .  (3) 

It  is  zdso  easy  to  show  that  the  error  propagation  operator  of  this  mul- 
tiplicative Schwarz  method  is 

{I  -  P2){I  -  Px), 

This  algorithm  can  therefore  be  viewed  as  a  simple  iterative  method  for 
solving 

with  an  appropriate  right-hand  side  5/,  . 

This  operator  is  a  polynomial  of  degree  two  and  thus  not  ideal  for  par- 
allel computing,  since  two  sequential  steps  are  involved.  If  more  than  two 
subspaces  are  used,  this  effect  is  further  pronounced,  even  if  the  degree  of 
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the  polynomial  representing  the  multiplicative  algorithm  often  is  lower  than 
maximal.  This  is  so  because  a  product  of  two  projections  associated  with 
subregions  which  do  not  overlap,  vanishes;  cf.  the  discussion  in  Widlund 
[31].  The  basic  idea  behind  the  additive  form  of  the  algorithm  is  to  work 
with  the  simplest  possible  polynomial  in  the  projections.  Therefore  the 
equation 

Puh  =  (Pi  +  P2  +  •-  +  PN)uh  =  g'h  ,  (4) 

is  solved  by  an  iterative  method.  Since  the  operator  P  is  symmetric,  with 
respect  to  the  bilinear  form,  and  positive  definite,  the  method  of  choice  is 
the  conjugate  gradient  method.  Equation  4  must  have  the  same  solution  as 
equation  1,  i.e.  the  correct  right-hand  side  must  be  found.  Since  by  equation 
1,  a{uh,<f>h)  =  f{4>h)  ,  the  right-hand  side  gj,  can  be  constructed  by  solving 
equation  3  for  all  values  of  i  and  adding  the  results.  It  is  similarly  possible  to 
apply  the  operator  P  of  equation  4  to  any  given  element  of  V^  by  applying 
each  projection  Pi  to  the  element  and  adding  the  results.  Most  of  the  work, 
in  particular  that  which  involves  the  individuaJ  projections,  can  be  carried 
out  in  parallel. 

We  now  describe  the  additive  Schwarz  method  introduced  in  Dryja  and 
Widlund  [14];  cf.  also  Dryja  [11]  .  We  start  with  the  same  triangular  (sim- 
plicial)  nonoverlapping  substructures  fi,  1  that  have  been  considered  before. 
Since  Schwarz-type  domain  decomposition  algorithms  use  overlapping  sub- 
regions,  we  extend  each  substructure  to  a  larger  region  fij  .  We  assume  that 
the  overlap  is  generous  assuming  that  the  distance  between  the  boundaries 
dVti  aJid  on',  is  bounded  from  below  by  a  fixed  fraction  of  H,,  the  diameter 
of  ili  .  We  also  assume  that  dVl'^  does  not  cut  through  any  element.  We 
make  the  same  construction  for  the  substructures  that  meet  the  boundary 
except  that  we  cut  off  the  part  of  Q[  that  is  outside  of  H  . 

The  analysis  of  Schwarz  methods  is  more  complicated  when  the  bound- 
aries of  the  different  subdomains  Q[  intersect  at  one  or  several  points;  cf. 
the  discussion  in  Lions  [20].  Such  a  situation  occurs  if  the  region  is  L-shaped 
and  is  partitioned  into  two  overlapping  rectangles.  We  discuss  such  a  more 
complicated  situation  in  Section  5  . 

Our  finite  element  space  is  represented  as  the  sum  of  N+1  subspaces 

V''  =  Vq^  +  V^^  +  ...  +  V^  . 

The  first  subspace  Vq  is  equal  to  V^ ,  the  same  coarse  global  space  of  contin- 
uous, piecewise  lineaj  functions  on  the  coarse  mesh  defined  by  the  substruc- 
tures ili  that  we  introduced  in  Section  2.  The  other  subspaces  are  related  to 
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the  subdomains,  in  the  same  way  as  in  a  traditional  Schwarz  algorithm,  i.e. 
VjA  _  V^  f\HQ{Q[).  The  computation  of  the  projection  of  an  arbitrary  func- 
tion onto  the  subspace  V^  involves  the  solution  of  a  standard  finite  element 
linear  system  of  algebraic  equations  which  is  on  the  order  of  N.  This  coarse, 
global  approximation  of  the  elliptic  equation  is  of  the  same  type  as  the  local 
problems  associated  with  subdomains.  The  only  real  difference  between  the 
problem  related  to  the  first  subspace  and  the  others  lies  in  the  way  that  the 
right  hand  side  of  the  linear  system  is  generated  as  weighted  averages  with 
weights  determined  by  the  basis  functions  associated  with  the  coarse  mesh. 
We  note  that  if  we  maJce  the  dimension  of  all  the  subspaces  approximately 
equal,  then  we  will  have  N  +  1  linear  systems,  each  with  about  than  N 
unknowns  ,  to  solve  in  each  step  of  the  iterative  solution  of  a  linear  system 
with  about  iV^  unknowns. 

It  is  well  known  that  the  number  of  steps  required  to  decrecise  an  appro- 
priate norm  of  the  error  of  a  conjugate  gradient  iteration  by  a  fixed  factor  is 
proportional  to  y/n  ,  where  k  is  the  condition  number  of  P  ;  see  e.g.  Golub 
and  Van  Loan  [16]  .  We  therefore  need  to  establish  that  the  operator  P 
of  equation  4  is  not  only  invertible  but  that  satisfactory  upper  and  lower 
bounds  on  its  eigenvalues  can  be  obtained.  A  constant  upper  bound  can  eas- 
ily be  obtained  for  P;  cf.  Section  4.  For  certain  other  additive  algorithms  a 
useful  technique  is  based  on  strengthened  Cauchy  inequalities;  cf.  Mandel 
and  McCormick  [21],  Widlund  [32]  and  Yserentant  [34]. 

A  lower  bound  can  often  conveniently  be  obtained  by  using  a  lemma, 
given  by  Lions  [20]  for  the  case  of  iV  =  2;  a  proof  is  also  given  in  Widlund 
[32]. 

Lemma  1  Let  u/,  =  X).=i  "A.n  where  Uh,i  €  Vj,  fee  a  representation  of  an 
element  of  V^  =  Vi  +  ...  +  Vjv.  //  the  representation  can  be  chosen  so  that 
EiLi  ainh,„Uh,i)  <  C^a{uH,Uh)yuH  G  V\  then  A„.„(P)  >  C^^ ■ 

We  remark  that  it  is  clear  that  Co  decreases  if  the  we  expand  the  sub- 
spaces.  This  follows  from  the  fact  that  there  is  a  larger  choice  in  selecting 
«h,i  G  Vf.U  we  can  expand  the  subspaces  without  worsing  the  upper  bound, 
and  that  is  often  possible,  our  estimate  of  k(P)  improves.  On  the  other  hand 
a  larger  subspace  also  means  that  the  subproblems  have  more  variables  and 
that  they  are  worse  conditioned.  For  the  special  case  of  the  classical  Schwarz 
method  on  two  regions,  this  tradeoff  is  well  understood;  for  a  discussion  of 
precise  estimates  of  the  rate  of  convergence  cf.  Bj0rstad  and  Widlund  [3]. 

As  was  previously  pointed  out,  the  framework  with  subspaces  and  projec- 
tions has  already  proven  quite  useful  not  only  for  the  study  of  Schwarz-type 
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methods,  but  also  for  iterative  refinement  methods.  Similarly,  it  can  be 
shown,  that  the  bound  for  the  condition  number  of  the  iteration  operator  of 
the  hierarchical  basis  multigrid  method,  introduced  by  Yserentant  [34]  can 
be  derived  by  using  these  techniques.  Yserentant's  algorithm  involves  the 
use  of  a  direct  sum  of  cleverly  chosen  subspaces,  which  are  quite  different 
from  those  considered  here.  All  except  one  of  the  subproblems  that  result 
are  very  well  conditioned.  However,  since  a  direct  sum  of  subspaces  is  used, 
there  is  no  flexibility  in  the  representation  of  the  elements  of  F''.  This  fact 
adds  to  the  understanding  why  Yserentant's  method  is  much  less  attractive 
in  the  case  of  three  dimensions  since  the  best  possible  Co  and  the  condition 
number  grows  rather  rapidly  in  that  case,  when  the  mesh  is  refined. 

4     Analysis  of  an  Additive  Schwarz  Method. 

In  this  section,  we  will  study  the  method  introduced  in  the  previous  section 
cind  give  a  proof  of  the  following  result. 

Theorem  1  The  operator  P  of  the  additive  algorithm  defined  by  the  spaces 
V^  and  V/^  satisfies  the  estimate  k{P)  <  const. 

Here  as  elsewhere  in  this  paper,  the  constants,  in  our  estimates,  are 
independent  of  h  and  H. 

An  upper  bound  for  the  spectrum  of  P  is  quite  easy  to  obtain.  We  only 
have  to  note  that  for  i  >  1  , 

a(P,Uh,Uh)  =  a{PiUh,P,Uh)  =  an»(P,u/i,  PiU/j)  <  aQ>{uh,Uh). 

We  recall  that  the  subscript  indicates  the  domain  of  integration  of  the  bi- 
linear form.  The  basic  observation  here  is  that  P,u;,  can  be  regarded  as  a 
projection  of  n^{n'^)f]V^  onto  n^(n[)f]V''  .  In  the  partion  of  the  region 
considered  here,  each  point  is  covered  by  subregions  Q[  a  finite  number  of 
times  .  A  constant  upper  bound  of  the  eigenvalues  of  P  is  therefore  obtained 
by  noting  that,  additionally,  the  norm  of  Pq  is  equal  to  one. 

The  lower  bound  is  obtained  by  using  Lemma  1.  We  partition  the  finite 
element  function  Uh.  as  follows.  We  first  choose  Uhfl  £  V^  .  By  using 
smoothing  and  interpolation,  cf.  e.g.  Strang  [26],  we  can  find  a  linear  map 
In  into  V^ ,  which  is  bounded  in  .^o(n)  and  which  satisfies 

||«;i  -  iHUhWiim  <  const. n\uh\fmn)  ■  (5) 
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Let  Wh  -  Uh-  in^h  ■  The  other  terms  in  the  representation  of  Uh  are  defined 
by  Uh,i  =  IhiO,Wh)  ,i=  1,---,N  .  Here  Ih  is  the  interpolation  operator  into 
the  space  V'^  and  the  9,  define  a  partition  of  unity  with  Oi  6  Co'i^i)  and 
^^,(i)  =  1  .  Because  of  the  relative  generous  overlap  of  the  subregions, 
introduced  in  Section  3,  these  functions  can  be  chosen  so  that  V^,  is  bounded 
by  const. /Hi.  By  using  the  linearity  of /^ ,  we  can  easily  show  that  we  obtain 
a  correct  partitioning  of  u/, .  In  order  to  estimate  the  semi-norm  of  Uh,i  ,  we 
work  on  one  element  K  at  a  time.  We  obtain 

\uhAh(K)  <  2|^u'/.lw>(/C)  +  ^\hi{e,  -  ^)u'/,)iH.(K)  • 

Here  Ti  is  the  average  value  of  fl,  over  K.  It  is  easy  to  see,  by  using  an  inverse 
inequality,  that 

\hiio. -f,)wh)\„^K)  <  const.  h-'\\h{{e,  -el)w^)\\L,^J,■^ . 

We  can  now  use  the  fact  that  on  K,  0,  differs  from  its  average  by  at  most 
const,  hf Hi  .  After  summing  over  all  elements  of  il'-  ,  we  arrive  at  the  in- 
equality 

We  now  sum  over  all  i  and  use  that  each  point  in  il  is  covered  only  a  fixed 
number  of  times.  We  then  obtain  a  uniform  bound  on  Cq,  and  conclude  the 
proof  of  Theorem  1,  by  estimating  the  two  terms  of 

by  |i'/»l//i(n)  •  The  bounds  follow  by  using  the  boundedness  of///  in  H^  and 
inequcdity  5,  respectively. 

5      Iterative  Substructuring  Methods. 

In  this  section,  we  show  how  we  can  obtain  an  iterative  substructuring 
method  by  using  the  framework  with  subspaces  and  projections  developed 
in  Section  3.  We  primarily  consider  problems  in  the  plane.  We  give  an  es- 
timate of  the  condition  number  of  the  operator  P  ,  which  corresponds  to  a 
specific  choice  of  subspaces.  This  proof  can  be  carried  out  without  using  an 
extension  theorem.  We  then  indicate  how  results  for  the  the  special  case  of 
two  substructures  can  be  used  to  derive  different  algorithms,  among  them 
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one  due  to  Bramble,  Pasciak  and  Schatz  [5].  We  also  briefly  discuss  the 
three  dimensional  case. 

We  assume  that  the  region  is  divided  into  substructures  as  in  Section  2. 
An  additive  Schwarz  method  is  introduced  by  using  the  coarse  space  V^ 
and  subspaces  V-'j,  which  are  related  to  an  edge  of  a  substructure.  Thus  let 
r,j  be  the  edge  which  is  common  to  two  adjacent  substructures  fl,  and  fij. 
Then,  1/.)  =  n^{ni,)r\V'',  where  Jl„  =  n.Ur.,Ufi,  • 

Compared  with  the  subspaces  used  in  the  previous  section,  we  use  less 
overlap  in  the  sense  that  only  the  elements  of  V"  differ  from  zero  at  the 
vertices  of  the  substructures.  This  is  reflected  in  a  poorer  bound  on  the 
condition  number. 

Theorem  2   The  operator  P  of  the  additive  algorithm  defined  by  the  spaces 

V^  and  V,5  satisfies  the  estimate  k(P)  <  const. {I  +  log{nih)f  . 

By  using  the  same  method  as  in  Section  4,  we  can  easily  show  that  the 
eigenvalues  of  P  <  4  .  In  the  proof  of  the  lower  bound  of  the  spectrum  of 
P,  we  use  Lemma  1  and  the  following  Lemma,  which  plays  an  important 
role  in  the  more  traditional  theory  for  iterative  substructuring  algorithms. 
Variations  of  this  result  ,  which  dates  back  at  least  to  1966  ,  are  given  in  a 
number  of  papers  ;  see  e.g.  Bramble  [4]  ,  Bramble  ,  Pasciak  and  Schatz  [5] 
or  Yserentant  [34]  . 

Lemma  2  Let  a  be  any  value  ofuh{x),  with  x  £  Hi.  Then 

ll"A  -  a|lS,oo(n.)  <  <^onst.  (1  +  log(.ff//i))|u/,|^,jQ.j  . 

We  note  that  this  result  holds  only  for  regions  in  two  dimensions  and 
that  it  resembles  a  Sobolev  inequality.  Since  all  the  elements  of  V^^  vanish 
at  the  vertices  of  the  substructures,  and  we  must  pick  the  interpolant  ///u;, 
as  the  element  in  V^  in  the  representation  of  u/,  .  It  is  easy  to  show  that 
|///U/,|?^,  Q     can  be  estimated  by 

Y,{uH{Vk)-UHm)\ 

where  the  Vj^s  are  the  vertices  of  the  substructure  fi^  .  It  then  follows  from 
Lemma  2  that 

l^//«/.l^j(n.)  <  const.{\  +  log(F//i))|u/,|^,jj|.j  .  (6) 

Let  «;/,  =  «/i  —  In^h-  As  in  the  previous  section,  we  use  a  partion  of 
unity.    The  argument  is  now  more  complicated  since  there  is  less  overlap. 
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The  elements  of  V,^,  used  in  the  representation  of  u/,,  are  given  by  the 
formula 

Uh,ij  =  Ih{6ijWh)  ■ 

Since  Wh  vanishes  at  the  vertices  of  the  substructures  and  we  only  use  values 
of  dij  at  nodal  points,  we  only  need  a  partion  of  unity  at  the  nodal  points 
which  are  not  substructure  vertices.  The  function  ^,_,  must  be  equal  to  1  at 
all  nodal  points  in  the  interior  of  T.j  ,  since  all  the  other  cut-off  functions 
vanish  there,  and  it  must  vanish  at  the  corresponding  nodes  on  the  other 
edges  of  fi,  and  fi^  .  It  is  easy  to  see  that  V^.j  therefore  must  grow  as  fast 
as  const,  /r  ,  where  r  is  the  distance  to  the  closest  endpoint  of  F.j.  Cut-off 
functions,  for  which  this  is  also  an  upper  bound,  can  indeed  be  constructed. 
The  necessary  bound  is  obtained,  one  substructure  at  a  time.  In  order 
to  estimate  |u/i,,jlHi(n.)  '°  terms  of  |u'A|//i(n.)'  ^^  ^^^^  consider  the  few 
elements  of  f2,  for  which  an  end  point  of  F^j  is  a  vertex.  A  straightforward 
calculation  shows  that  the  contributions  to  |"/i,<j|^i(n,)  from  these  triangles 
are  smaller  than  those  to  |«^/i|//i(n  \-  We  now  use  arguments  similar  to  those 
of  the  previous  section  to  obtain  an  estimate  of  the  integral  over  the  rest  of 
the  substructure.  Working  with  one  element  at  a  time,  we  obtain 

By  using  the  fact  that  Wh  does  not  change  if  a  constant  is  added  to  u^,  the 
formula  for  u/i,,j  and  that  |lu>/i||£,^(;c)  <  2||u;,||^^,j^  >   we  obtain 

Wh,ij\Hi(K)  ^  '^l^hlnnK)  +  const.{h/r)^\\uh  -  «lli^(fi,)  • 

The  sum  of  the  first  term  over  the  elements  can  be  estimated  by  using 
inequality  6.  Since  the  number  of  elements  decreases  to  zero  linearly  with 
the  distance  r,  the  sum  of  the  second  expression  over  all  the  elements  of  the 
substructure,  except  those  at  the  endpoints  of  F.j,  caji  be  estimated  by 


rH 


Here  H  represents  the  diameter  of  fi,  and  h  the  minimum  distance  of  any 
other  nodal  point  to  the  end  points  of  F,j  .  The  estimate  on  Cg  and  the 
proof  of  the  whole  theorem  is  now  concluded  by  using  Lemma  2. 

This  algorithm,  introduced  and  analyzed  as  a  method  based  on  subspaces 
and  projections,  can  equally  well  be  understood  in  terms  of  a  splitting.  Let 
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us  consider  the  case  where  the  components  of  the  right  hand  side  of  the 
system  of  equations  corresponding  to  the  interiors  have  been  set  to  zero  in  a 
preliminary  step.  The  right  hand  side  of  a  linear  system  related  to  n,_,  and 
V/^  then  differs  from  zero  ordy  on  T.j.  The  Schur  complement  associated  with 
this  edge  can  be  shown,  straightforwardly,  to  be  the  sum  of  the  corresponding 
blocks  in  the  four  by  four  block  representation  of  the  Schur  complements  5^'' 
and  S^^\  which  were  introduced  in  Section  2  . 

One  of  the  attractive  features  of  the  framework  first  introduced  in  Section 
3  is  the  ease  by  which  the  subproblems  can  be  replaced  by  preconditioners. 
Let  us  consider  the  additive  Schwarz  method  for  the  problem  discussed  in 
the  beginning  of  Section  3.  We  can  write  the  projection  Pi  in  matrix  terms. 
After  a  suitable  permutation  of  the  variables,  it  is  seen  to  correspond  to 

It  is  easy  to  see  that  this  matrix  is  symmetric  in  the  K-  inner  product 
which  corresponds  to  the  bilinear  form.  K  K^^^  ,  and  the  other  matrices 
which  play  similar  roles,  are  replaced  by  inverses  of  preconditioners  for  the 
subproblems,  then  it  is  easy  to  see  that  the  resulting  algorithm  converges  and 
that  its  condition  number  can  be  estimated  immediately  in  terms  of  k(P) 
and  bounds  for  the  local  preconditioners.  Using  this  method,  a  number 
of  algorithms  can  be  derived  from  a  basic  method  based  on  projections 
and  subspaces.  Thus,  if  the  Schur  complement  corresponding  to  problem 
defined  on  ftij  is  replaced  by  the  square  root  of  the  discrete,  one-dimensional 
Laplacian,  denoted  by  Iq^^^  in  Bramble,  Pasciak  and  Schatz  [5],  we  obtain 
the  main  algorithm  of  that  paper.  The  estimate  of  the  condition  number 
of  the  resulting  method  can  be  obtained,  by  using  the  argument  just  given, 
combining  Theorem  2  and  a  bound  for  the  condition  number  of  problems 
defined  on  the  union  of  two  subregions. 

We  conclude  this  paper  by  a  remark  on  the  part  of  preconditioner  which 
corresponds  to  the  coarse  global  problem.  In  two  dimensions,  we  have  used 
the  subspace  V^  for  this  purpose.  The  related  quadratic  form  has,  in  a 
special  case,  the  form 

It  is  easy  to  show  that  this  form  equally  well  could  be  written  as  a  double 
sum  over  (u/i(Vfc    )  — «),  )^,  where  u]^'  is  the  average  value  of  the  three  values 
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associated  with  the  vertices  of  fi,.  In  three  dimensions,  the  corresponding 
quadratic  form,  with  the  sums  and  averages  calculated  with  respect  to  the 
variables  associated  with  the  wire  baskets  of  the  individual  substructures, 
provide  an  important,  global  part  of  the  powerful  iterative  substructuring 
methods,  which  have  been  developed  by  Bramble,  Pasciak  and  Schatz  [8] 
and  Dryja  [12]  . 
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